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ABSTRACT 


Two primitive variable, pressure based, flux-split, RNS/NS solution procedures for 
viscous flows are presented. Both methods are uniformly valid across the full Mach 
number range, Le., from the incompressible limit to high supersonic speeds. The first 
method is an "optimized" version of a previously developed global pressure relaxation 
RNS procedure. Considerable reduction in the number of relatively expensive matrix 
inversion, and thereby in the computational time, has been achieved with this procedure. 
CPU times are reduced by a factor of 15 for predominantly elliptic flows (incompressible 
and low subsonic). The second method is a time-marching, 'linearized' convection 
RNS/NS procedure. The key to the efficiency of this procedure is the reduction to a single 
LU inversion at the inflow cross-plane. The remainder of the algorithm simply requires 
back-substitution with this LU and the corresponding residual vector at any cross-plane 
location. This method is not time-consistent, but has a convective- type CFL stability 
limitation. Both formulations are robust and provide accurate solutions for a variety of 
internal viscous flows to be provided herein. 
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1. INTRODUCTION 


Various asymptotic approximations to the complete Navier-Stokes (NS) equations 
have been used to provide detailed, efficient and accurate flowfield descriptions for a 
significant class of large Reynolds number (Re) flows [1-5]. If these approximations are 
combined in a single system of composite equations, the resulting time-dependent system 
is termed the Reduced Navier-Stokes approximation (RNS). The RNS system in 
appropriate (£,T|,Q coordinates is such that only streamwise or t, diffusion terms in the full 
NS equations are higher order and therefore neglected throughout the flow domain. These 
terms are retained in a deferred corrector (DC) which may be recovered when necessary. 
The lowest order RNS approximation consists of the full Euler equations plus all the 
boundary-layer diffusion terms required to satisfy appropriate no-slip boundary conditions 
on various solid boundaries. In this manner, all acoustic (elliptic) influences contained in 
the full Navier Stokes (NS) equations, are retained. The resulting RNS system allows tor 
upstream or elliptic influence and contains all the dominant physics associated with large 
Reynolds number strong viscous-inviscid interactions. A pressure based flux-splitting 
procedure is applied to the Euler component ot these equations. This leads to a global 
relaxation procedure for the pressure, and for velocities in reverse flow regions. The 
convective and acoustic fluxes are treated independently. Theretore, the appropriate 
domain of dependence is automatically represented by the differencing of the convective 
and acoustic (pressure) gradients. 

Two efficient, primitive variable RNS/NS solution procedures, valid across the full 
Mach number range, are presented in this paper. The primary features of these solvers are 
(a) that the boundary conditions and discretization procedure are controlled by the physics 
of the of the lowest order terms of the RNS system; (b) that the higher order diffusion. DC 
terms, in the NS system can be introduced if necessary in the final RNS inversion to obyain 
the full NS solution; (c) that they possess sharp shock capturing properties; within three 
grid points; (d) unlike other Euler-based Navier-Stokes methods, that become 
unconditionally unstable in the incompressible limit, these formulations pertorm efficiently 
for all Mach numbers from incompressible to high supersonic and; (e) a pressure velocity 
flux-split discretizadon is implemented in both of these solvers. 



The first of the RNS/NS solvers is an optimized version of a global pressure 
relaxation RNS procedure previously presented [1-3]. The main drawback ot the original 
algorithm is the cumulatively high number of expensive LU of cross-plane coefficient 
matrix inversions that are required. At least three non-linear Newton iterations or LU 
inverions are needed at each cross-plane. For three-dimensional (3D), turbulent flows, 
wherein fine grids are needed to resolve the thin boundary layers, the matrix inversions 
become computationally expensive and impacts adversely on the efficiency and speed of 
the solution procedure. One of the primary goals of the present study is to minimize this 
cost. Toward this end various optimizational techniques are adopted: (a) recast the 
governing equations in "delta" form, so that a single LU inverse at a given cross-plane can 
be applied as an approximate LU for several subsequent stations downstream; (b) global 
under-relaxation of flow variables, with time terms in the governing equations (c) local 
under-relaxation of the flow variables during the non-linear iterations and; (d) a judicious 
initialization of flow variables at a given cross plane with those computed at a preceding 
station in the prior sweep. With these techniques considerable reduction in the required 
number of LU inversions, and consequently in the computational time results. For 
instance, the optimized version of the code runs about 15 times faster than the original 

"unoptimized" procedure for laminar, incompressible flow in a 90° curved duct, and a 
33x41x41 grid. The optimized code requires a single LU inverse to attain a steady-state 
solution. This compares with 174 inverions per sweep for the original version. In this 
particular case. 32 global sweeps are required to converge the maximum residuals to the 
prescribed tolerance level. Note, that it requires exactly "imax-1" number of global sweeps 
for incompressible flows to converge to machine accuracy; where "imax" is the number of 
stations in the streamwise direction, as the downstream pressure boundary condition 
traverses a single grid cell for even,' pass. 

The second algorithm, i.e., linearized 1 convection model, is a time-marching RNS 
procedure. This is a straightforward algorithm that involves just a single LU inverse at the 
inflow cross-plane. The method marches in pseudo-time with the application of a single 
back-substitution, for the LU inverse and the corresponding residual vector at each cross- 
plane. during a given sweep. With this mathematical operation, all the flow variables are 



marched a single time-step in a given sweep. Note, that this procedure is not time- 
consistent. This process of time-marching, is repeated until a steady-state is attained. A 
convection-only CFL type time-limitation, discussed in reference [4-5] for the time 
consistent algorithm, is applicable for this method. Since flow variables at all cross-plane 
locations can be marched a single time-step simultaneously i.e., they are fully uncoupled 
numerically in the axial or flow direction, this procedure lends itself to parallelization. This 
method is particularly efficient for incompressible and low subsonic Mach numbers. For 
supersonic flows, the "optimized 1 pressure relaxation procedure is more suitable and 
therefore preferred. 

Both methods implement a pressure velocity flux-split discretization [1-3]. The 
RNS/NS system of equations is quasi-linearized and discretized as described previously in 
[4,5]. A system of simultaneous algebraic equations (for the cross-plane in three 
dimensions) is solved with a sparse matrix direct solver (SMDS). The use of this solver is 
dictated by robustness and consistency consideradons discussed in previous papers [1-3]. 

Both methods have been validated through a series of three-dimensional internal 
flow configurations: (a) Laminar and turbulent, incompressible flow in an S-shaped duct; 
(b) laminar and turbulent, compressible flow in a symmetric, square cross-section choked 
convergent-divergent nozzle and, (c) turbulent, supersonic flow in a generic inlet 
geometry. 

The algebraic Baldwin-Lomax eddy viscosity model, modified to include multiple 
wall effects, is employed for turbulence closure. 

2. GOVERNING EQUATIONS 

The reduced form of the Navier-Stokes (RNS) equations is employed for both of 
the flow solvers considered herein. All higher order axial diffusion terms are retained only 
in the deferred-corrector (DC). The DC is explicitly introduced in the final RNS inversion 
to obtain a full NS solution, see ref. [10] 

In order to obtain the final RNS system of equations, the full NS equations are 
transformed from Cartesian into non-orthogonal generalized curvilinear coordinates in 
strong conservation form. These are as follows: 
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Y-momentum equation: 
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Z-momentum equation: 
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Energy equation ( constant total enthalpy): 
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Equation of State: 
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By eliminating temperature T from the equations (5) and (6), the following relationship is 
obtained: 
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It is this form of combined energy equation and the equation of state that is used in the 
present study. In the incompressible limit, the non-dimensional density is set to umty. The 
viscous terms, viz, r;’, r’ , T? , tj , etc., appearing on the right-hand-side (RHS) in 

equations 1-4 have been discussed previously in [6]. The definitions of the metric 
quantities, e.g., q x . ?v, tix, Cv” etc ” ^ the contravariant velocities U . V. and vv 

appearing in these equations have also been defined in [6]. The RNS set of equations are 
obtained from the full NS equations by simply dropping the higher order diffusion (q- 
derivatives) terms. These are retained in the DC. The final form of the momentum 
equations used in these solvers are obtained by taking the covariant momentum balances in 
the %,r|,C directions. 



3. DISCRETIZATION 


The pressure velocity flux-split technique, originally developed by Rubin and Lin 
[1,7], is applied to the RNS equations in generalized non-orthogonal curvilinear 
coordinates. Fig. 1 depicts the discretization location of the equations and their appropriate 
groupings, at a typical cross-plane, of a given axial station The ^-momentum and the 

energy equations are discretized at the grid points. The continuity equation is cell centered 
and the r| and ^-momentum equations are located at the half points. All ^-derivatives in 
the continuity equation are backward differenced and the cross-flow (r|,Q derivatives are 
two-point trapezoidal differenced. All axial convection terms in the momentum equations 
are upwind or flux-vector differenced and the cross-flow convection terms are either 3- 
point central differenced or two-point trapezoidal differenced. This depends on the 
location at which these derivatives are discretized. All diffusion terms are second order 
accurate three-point central differenced. The streamwise pressure gradient in the q- 
momentum equation is flux split into positive and negative contributions in accordance 
with the pressure flux- vector splitting technique [1]: 
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where is located at the same grid location as the velocity u , , and is given, to second 


order in Aq , by 
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For constant stagnation enthalpy, the variable oi in the above equations is given as 
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boundary value or the elliptic component of the pressure gradient. The discretization of 
implies that the "unknown" pressure P^ at the axial (marching) location i, is staggered at a 


distance (l-co)A^ upstream of the velocity up vp and wp The pressure P i at the grid point i 
is given by eqn. (9). In the incompressible limit, o)=0, and the pressure at the grid point i is 
P i = P i + 1' For supersonic flows, co=l and P i = P i . By neglecting the negative flux 

contribution of the streamwise pressure gradient, Le. by retaining only the hyperbolic 
component in P^ 
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, the multiple-sweep global pressure relaxation RNS formulation 


reduces to the single-sweep PNS (Parabolized Navier-Stokes) formulation. The 
discretization of cross-flow pressure gradients P^ and P^ in the r|-momentum and the 2- 
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momentum equations is quite straightforward and is given by two-point trapezoidal 
differencing. 



Fig.l Implementation of boundary condition in a q, cross-plane for an internal flow 

geometry. 


4. BOUNDARY CONDITIONS 


The appropriate boundary conditions, consistent with the discretized system of RNS 
equations, are as follows: 

inflow boundary \ 'z - g m ) : 

In the "streamwise" or £ direction, at the inflow boundary, all flow variables or 
derivatives are prescribed as follows: 


ytj 



« = u in (rj,C) 

v = v in^O 
w = w in (rJ ’^) 

M = Mj (zero for incompressible flow) 

P = P in (r 7’0 

= 0 (for internal, compressible flow) 

For incompressible flow calculations, the inflow pressure is not prescribed, but is 
calculated from Rj = P t . For subsonic flows, the pressure j as defined by eqn. (9), 

still retains some influence of the prescribed inflow solution Pj./. In these calculations the 
value of P { _ j changes to reflect the upstream pressure influence, see Rosenbaum and 
Rubin [8]. This allows for adjustment of inflow mass through changes in Pj and hence p( n . 

outflow boundary (q = q oul ) : 

At a subsonic or incompressible outflow, 0 < M * < 1 , without flow reversal, only the 

pressure or pressure gradient is required. For the PNS "supersonic" step, the negative 
pressure fluxes are not present, and therefore, the pressure is calculated at this boundary. 
Due to the presence of a "subsonic" boundary layer region, even in supersonic viscous 
flow cases, the pressure at the exit boundary, computed from the PNS step, is prescribed 
as the outflow pressure boundary condition for subsequent RNS steps. 

For zero flow reversal at the outflow, the velocities are calculated from the solver. 

For small reverse flow at the exit, the velocities can be prescribed from available 
experimental data, or computed by neglecting the negative velocity flux at the outflow 
boundary. This is similar to a FLARE approximation used in boundary layer theory. For 
the present study, the FLARE approximation is applied only at the outflow boundary. 



wall boundary: 

No slip and zero-injection (u=0, v=0, w=0) for viscous How computation. 

Zero-vorticity condition for inviscid/Euler computation. 

Pressure is computed at all wall grid points using the special regroupings discussed 
in [4,5], 

The density is computed at all wall points from the constant stagnation enthalpy 
condition. 

free boundaries ( external flow): 

The inplane velocities, pressure and density are specified. The velocities normal to the free 
boundary are computed. 

5. SOLUTION PROCEDURE 

For both procedures considered herein, the discrete system of fully coupled RNS 
system in generalized, non-orthogonal, curvilinear coordinates, at a given cross-plane, is 

written in the following delta form: = (r) rt ; where [a] represents the sparse 

cross-plane coefficient matrix; represents the solution vector in delta form and 

( r) n is the right-hand-side residual vector, evaluated at a previous time or iteration. Since 
the primary goal is to drive the local cross-plane residuals to a prescribed tolerance level, 
an "exact" inverse of W» not required. An approximate LU that is "close-enough ' in 

character to the exact one is found sufficient. In other words, in delta form, an exact LL 
of the coefficient matrix at an arbitrary cross-plane, can successfully drive the residuals, at 
several subsequent "streamwise" stations, to an acceptable level of accuracy. Furthermore, 
if the "optimizational" techniques, mentioned in the Introduction ot this paper, are 
implemented appropriately, then just a single LU ot the coefficient matrix, at the inflow 
cross-plane, has been found to be sufficient for the entire computation. Although, this 
misht result in an increased number of back-substitutions to resolve local non-linearities, 
the overall cost is considerably reduced. Back-substitution operations, with an 



approximate LU, is applied iteratively until the local non-linearities are resolved and the 
local residuals have achieved a prescribed tolerance level. If during this iterative back- 
substitution process the residuals show signs of growth, or if the local convergence ot 
non-linearities become excessively slow, a new LU is initiated. In addition, the number ot 
required LUs can be further reduced with the aid of under-relaxation, both globally 
(inclusion of time-terms) and locally. This process of back-substitution and LU 
decomposition is continually repeated until there is local convergence of all residuals at all 
the axial stations, in a given sweep. Multiple sweeps are generally required to converge all 
residuals globally. 

For the linearized convective time-marching RNS procedure, only one LU inverse 
at the inflow cross-plane is required. This linearizes the non-linear terms about the 
prescribed inflow conditions. Unlike the global pressure relaxation procedure, the non- 
linearities are not resolved locally but in a global fashion. A single back-substitution 
operation is required at any given axial cross-plane for all axial sweeps. As a result, all the 
flow variables in the computational domain evolve in pseudo-time through a single time- 
step. This process is repeated until steady-state conditions are attained. A convection-onlv 
CFL limitation is applicable, see ref. [4-5]. 

Fig. 2 depicts the comparison of convergence histories for the Unearned convection, time- 
marching model and the fully- imp Ucit, pressure relaxation procedure. The flow' problem 
under consideration is an incompressible, laminar flow in an S-duct at Re=790. The grid 
has 31 stations in the axial direction and 21x21 in a cross-plane. As seen in the figure, the 
convergence history for the linearized convection model shows a more gradual downward 
trend than the pressure relaxation model. However, the computational time required by 
both methods to converge to a steady-state is about the same i.e., approximately 50 
minutes on a RISC/6000 workstation. As noted earUer, apart from the difference in time- 
steps used in the two methods, the linearized convection model does not converge non- 
linearities at each axial location. This is required in the pressure relaxation procedure. The 
gain in time by the former method is clearly outweighed by the need for a larger number of 
axial sweeps or "time-steps" for convergence. 



Comparison of convergence histones for Uneanzed convection & pressure relaxation RNS procedures 



Fig. 2 Comparison of convergence histories for the Linearized convection model and 
the pressure relaxation procedure for a laminar, incompressible flow in an S-duct. 


Both procedures are quite robust and have been implemented in the same code. 
This provides the flexibility to switch from one algorithm to the other for a given flow 
problem in a single run. 


6. RESULTS 

In order to validate these RNS/NS flow solvers, a variety of 3D internal flow 
configurations, ranging in speed from the incompressible limit to high supersonic Mach 
number, have been investigated. An algebraic, zero-equation Baldwin-Lomax model, 
modified to handle multiple walls, has been applied in all the turbulent flow cases. 

Incompressible flow in an S-shaped duct: Incompressible flow in an S-shaped duct of 
constant area square cross-section is investigated. Two flow Reynolds numbers are 
considered, Re=790 (laminar) and Re=40,000 (turbulent). A combination of the linearized 
convection and the 'optimized' global pressure relaxation solvers is applied. The numerical 
results are compared with the measurements of Taylor et al. [9], The secondary flow 
phenomena in a S-duct is mainly pressure driven because of the very smooth-bend in the 



walls. Axial How separation is not observed. Accurate prediction of the boundary layer, 
even in the laminar case, is required to capture the small but complex secondary Hows. 
The grid in the cross-plane must be adequately refined in order to correctly predict this 
behavior. Grids with 33 axial stations and with 31x31 or 41x41 points in the cross-plane 
are prescribed for the laminar and turbulent cases, respectively. The symmetry plane grid 
for the turbulent case is shown in Fig. 3 The linearized convection code is applied in the 
initial phase of the computation to drive the maximum residuals down to 0.01. Thereafter, 
the "optimized" fully implicit relaxation solver, with infinite time step, is automatically 
initiated. In this mode, the residuals reach the prescribed levels of tolerance, typically 
0.0001, relatively fast. The convergence history of the laminar computation is depicted in 
Fig-4. 


Turbulent, Incompressible S-duct (Re=40,000) 
Grid (33x41x41) along the mid-plane 



Figure 3. Typical view of the geometry and the grid in the symmetry plane 



Convergence history for an incompressible S-duct 



Fig 4. Convergence history for a laminar, incompressible flow in an S-duct 


The streamwise velocity profiles, for both laminar and turbulent Hows, on the 
symmetry/midplane of the S-duct, is compared with the measured data of Taylor et al. [9]. 
Results at several axial stations, shown in Figs. 5(a-e) and Figs. 6(a-e), are depicted. There 
is an excellent agreement with the experimental data for both laminar and turbulent flows. 
This computation requires about 4 hours for the laminar and about 6 hours for the 
turbulent case on the RISC/6000 machine. This is roughly equivalent (without 
vectorization) to 30 and 50 minutes, respectively, on a Cray YMP supercomputer. 

In order to demonstrate the effect of the reduction in the required number of LU 
inversions on the computational time, a laminar incompressible case was considered. A 
single LU inversion is required by the "optimized" version of the global pressure relaxation 
solver, while 174 LUs are required by the original code during each axial sweep. Since the 
number of global sweeps approximately equals the number ot nodes in the axial direction, 
this results in an enhancement of almost 15. 



Figures 5(a-e). Comparison of the laminar axial velocity profiles on the symmetry 

plane at stations 1-5 along an S-duct 
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Figures 6(a-e). Comparison of the computed and measured turbulent axial velocity 

profiles for an S-duct 



(b) Turbulent , 3D convergent-divergent choked nozzle: Turbulent flow in a symmetric, 
convergent-divergent nozzle, with square cross-section, and an area ratio. 
Ainlet/Athroat=3-41, is computed with the linearized convection procedure. The Reynolds 

number is Re=5xl(P based on the throat diameter. This computation is tor design point 
conditions, so that the pressures are prescribed at both ends of the nozzle in accordance 
with the results of quasi- ID analysis. A Mach number of 0.2 is prescribed at the inlet. As 
expected, the nozzle is choked with a Mach number of unity at the geometnc throat. At 
the exit cross-plane, an average Mach number of 2.8 is computed. This is in excellent 
agreement with quasi- ID thoery. A grid with 31 axial stations and a highly stretched 
41x41 mesh in the cross-plane is employed. Figs. 7-8 depict the Mach number and pressure 
contours on a symmetry plane. The computed Mach number ratio between the inlet and 
the exit are in close agreement with quasi- ID analysis. 



Figure 7. Mach number contours on a symmetry plane 



Figure 8. Non-dimensional pressure contours on a symmetry plane 







( c ) Effect of back-pressure on a laminar, compressible flow through a choked 
convergent-divergent Nozzle: In order to study the effect of the back-pressure, 'p b \ on 

the performance of a convergent-divergent nozzle, a laminar flow at Re=500 is 
considered. The nozzle is identical to the one considered for the turbulent case. The grid 
has 51 axial stations and 33x33 nodes in a cross-plane. Mild stretching, of about 1.1 at the 
walls, is applied. This was adequate to resolve a relatively thick boundary layer associated 
with the lower Re. In this study, a series of computations with different back pressures at 
the exit was carried out. The first calculation was for a steady-state solution at design 
conditions, i.e., Pb=p e (complete expansion). Here 'p e ' denotes the exit pressure that 

would be attained under design conditions. From quasi-one-dimensional theory, the 
pressure ratio Pe/p^ for design conditions is approximately 27 for an area ratio 

Ae^throat ot: 3-51. Subsequently, the back-pressure is raised in a gradual fashion. First, 
the back-pressure is increased to p|,=3.6p e and a steady-state solution is obtained for this 

condition with the design-point solution as an initial guess. Next, the back-pressure is 
further increased to pb=6.25p e and again a steady-state solution is obtained with the 
previously obtained steady-state solution at p^=3.6p e applied as the starting solution. For 

this back-pressure, axial flow separation associated with oblique shock formation is 
predicted. In order to assess the effect of area ratio Ag/A^o^, the fl are angle of all the 

walls in the divergent section of the nozzles is then reduced from 5 to 3 degrees. This 
decreases the area ratio A e /A t | iroat from 3.51 to 2.3. The convergent section of the 

nozzle is not altered. The design point pressure ratio, Pe/pin, for the modified nozzle 
geometry is reduced to about 13.5 from quasi-one-dimensional analysis. The final imposed 
back-pressure, previously 6.25 times the design exit pressure in the original nozzle 
geometry, now reduces to 3.37 of the design exit pressure Le.. p{,=3.37p e . A steady-state 

solution is now obtained for the modified nozzle. For this calculation, the steadv-state 
solution associated with the original nozzle and pb=6.25p e . is applied as the starting 
solution. As expected the reverse flow regions, that were present in the original design no 



longer exists for the modified nozzle case. This provides a rapid tool for nozzle design 
studies. This computation required only 2 hours on RISC6000. 

Subsequent to the design point computation, only divergent portion of the nozzle 
need to be evaluated. As long as the nozzle remains choked any increment in the back- 
pressure does not affect the flow conditions upstream of the throat. Therefore, each 
steady-state solution, at higher back-pressures, was obtained relatively inexpensively. It 
took approximately 7 hours on the IBM RISC/6000 machine for the design condition 
computation. Each subsequent run took only about 2.5-3 hours. This translates to 
approximately 50 and 18-23 minutes, respectively, on a Cray YMP supercomputer. 

The effect of higher back-pressure on the pressure, skin friction, axial velocity and 
cross-flow velocity is depicted in Figs. 11-14. For Pb=Pe or Pe/Pin = 27, complete or 
optimal expansion takes place. For pb=3.57p e or Pe/pin = ^, A° w separation is 

imminent. This is visible from the skin friction coefficient along the lower wall in Fig. 12. 
Only a slight adverse pressure gradient is evident. As the back-pressure is further 
increased to pb=6.25p e or Pe/pin^, two strong oblique shocks are generated. Behind 

these shock waves the flow separates and the pressure increases to the imposed back- 
pressure. Since the core the flow is supersonic, the effect of increased back-pressure is 
propagated upstream through the subsonic boundary layers close to the wall. The skin 
friction and axial velocity clearly show the existence of the reverse flow region. Fig. 15 
shows the streamlines in the (x-y) symmetry plane, wherein reverse flow regions are 
clearly evident. Fig. 16 depicts the streamlines in the (x-y) symmetry plane for the modified 
nozzle. As expected, the shock induced reverse flow regions disappear. 



Effect of bac*.-tx«suf« (pj and divergence angle on the flow 



Figure 11 Comparison of pressure on the lower wall of the symmetry plane for 

various imposed back-pressures 


Effect of bacH-pre&surt (pj and divergence angle on die flow 



Figure 12 Comparison of skin friction coefficient Cf on the lower wall of the 
symmetry plane for various imposed back-pressures 
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Figure 13 Comparison of axial velocity "u" at the exit cross-plane for various 

imposed back-pressures 
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Figure 14 Comparison of cross-flow velocity "v" at the exit cross-plane for various 

imposed back-pressures 
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Figure 15 Streamlines on the (x-y) symmetry plane for the imposed back-pressure of 

Pb=6.25p e & Ag/Athroat of 351 
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Figure 16 Streamlines on the (x-y) symmetry plane for the imposed back- 
pressure of pb=3.37p e & Ag/Athroat of 2.3 

Grid Convergence and Accuracy: In order to assess grid convergence and accuracy of 
these calculations, three grids, 31x17x17, 31x33x33 and 51x33x33, with refinement in all 
three directions, were considered. Calculations were performed at design conditions for a 
laminar, compressible flow. Although Figures 17-20 suggest further grid refinement in all 
three directions to obtain grid independence, the solution obtained on the finest grid 
51x33x33 can be considered reasonably accurate. All plots, except pressure, show 
significant disparity between the solution obtained with a coarse 31x17x17 grid and that 
obtained with finer 31x33x33 and 51x33x33 grids. The pressure, however, is not affected 
by grid refinement. This is due to a predominantly inviscid nature of the flow. However, 
the accuracy of velocities and skin friction calculation is enhanced by clustering more grid 
points in the boundary layer. The effect of axial refinement is minimal on skin friction. It. 
however, increases the peak slightly of both, axial and cross-flow components of the 
velocity. The skin friction appears to be affected solely by refinement in the cross-plane. 

.After an initial decline in the skin friction, due to a developing boundary layer, it 
continuously increases throughtout the length of the nozzle. Both, skin friction and 
pressure plots, suggest an accelerating flow in the nozzle. Faster accleration rates are seen 
in the divergent section. A sudden expansion near the throat is evident in the pressure plot. 
This is also reflected in the increment of skin friction in the proximity of the throat. The 
axial velocity maximum occurs closer to the walls than the center line. This is seen in the 



axial velocity plot of Fig. 19. The reason being the existence of extremely thin boundary 
layers on the four adjacent walls, caused by a tremendous flow acceleration, that do not 
interact with each other in the central core region. As a result, a large portion of the core 
flow remains inviscid. 



Figure 17 Skin friction coefficient Cf on the lower wall of the (x-y) symmetry plane 
for 31x17x17, 31x33x33 and 51x33x33 grids 
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Figure 18 Non-dimensional pressure on the lower wall of the (x-y) symmetry plane 
for 31x17x17, 31x33x33 and 51x33x33 grids 




Figure 19 Axial velocity "u" at the throat for 31x17x17, 31x33x33 and 51x33x33 

grids 
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Figure 20 Cross-flow velocity "v" at the throat for 31x17x17, 31x33x33 and 

51x33x33 grids 


(d) Turbulent, Supersonic flow in a Generic Three-dimensional, 10° single-ramp inlet: 

Turbulent, supersonic flow, in a generic three-dimensional, 10° single-ramp inlet geometry 
is investigated with the 'optimized' solver. The inlet Mach number is M^^.O and the 

Reynolds number is Re=5xl0-- Efficient simulation of complex shock patterns in such 
internal configurations is examined. The grid in a symmetry plane is shown in Fig.2l. In 
order to study grid dependency and accuracy of the solution, computations were carried 
out on three different grids. .Although, grids I and II have same dimensions, i.e., 41 nodes 
in the axial direction and 41x41 in the cross-flow directions, the first node in the straight 



duct portion of the inlet for grid I is 0.002 away from the walls as opposed to 0.0008 in 
grid II. Grid III has the same dimensions and distribution of mesh points in the grid II 
cross-plane, but has 61 nodes in the axial direction. On average, the computation on a 
41x41x41 grid takes about 45 minutes per global pressure relaxation sweep on a 
RISC/6000 workstation. This is equivalent to about 6 minutes on a Cray YMP 
supercomputer. A total of about 6 global sweeps are required to converge all the residuals 
to the tolerance level of 0.0001. 
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Figure 21. Typical view of the grid in the symmetry plane 

Figs.22(a-d) depict the axial "u" and cross-flow velocity "v" profiles along the 
symmetry plane at axial stations 1 and 2. Solutions on all three grids have been used for 
this comparison. There are no noticeable changes in "u ", and for "v" there are only minor 
differences. 
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Figure 22(a-b). 


and "v" velocity profiles at "station 1 " along the mid-plane 



Figure 22(c-d). "u " and "v " velocity profiles at ", station 2 " along the mid-plane 

Figs.23(a-b) depict the non-dimensional pressure on the symmetry plane lower and 
upper walls. With reference to Fig. 23(a), the higher pressure level, behind the strong 
oblique shock, remains practically constant to the point where the ramp becomes 
horizontal (x=1.0). A steep decline in the pressure profile at this point is due to an 
expansion fan generated at the convex comer, x=1.0. Further downstream of this point, 
the pressure profile shows a slight tendency to rise. This is due to the impingement of the 
reflected oblique shock from the upper to the lower wall in the vicinity, x=1.9. Since the 
intensity ol the doubly reflected shock is considerably diffused at this point, the rise in 
pressure is comparadvely small. 



Figure 23(a-b). Non-dimensional pressure on the lower and upper walls, respectively , 

along the symmetry plane 



Fig. 23(b) depicts the pressure on the upper wall along the symmetry plane. In the 
region, 0.0 < x < 1.0, the pressure on the upper wall remains fairly constant. The strong 
oblique shock, that originates at the beginning of the ramp (x=0.0 at lower wall), impinges 
the upper wall at x=1.0. This causes an adverse pressure gradient in the region, 1.0 < x < 
1.4. If this negative pressure gradient is sufficiently high, it can cause flow reversal. 
Subsequently, in the range, 1.4 < x < 2.0 (exit), a steady drop in pressure is observed. This 
can be attributed to the influence of the expansion fan, that originates from x=1.0 on the 
lower wall and interacts with the reflected shock. This complex interaction of shock wave, 
expansion fan and boundary layer is seen in Figs.24(a-b). 
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Figure 24(a-b). Mach number and non-dimensional pressure contours on the 

symmetry plane 


7. CONCLUSIONS 

The validity of two primitive variable RNS/NS flow solvers, for calculations from 
the incompressible limit to supersonic speeds has been demonstrated. The required number 
of LU inversions has been significantly reduced for the global pressure relaxation 
procedure. As a result, a speed-up of almost 15 has been achieved for an incompressible 
laminar flow in a S-shaped duct. For non-separated supersonic flows, this method is even 
more efficient. The efficiency and the validity of the time-marching, linearized convection 
model presented herein, has also been validated with S-duct, and convergent-divergent 
choked nozzle computations. Agreement with the data has been quite good thoughout. 
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Implementation of boundary condition in a cross-plane for an internal flow 

geometry. 


Comparison of convergence histones for Linearized convection & pressure relaxation RNS procedures 



Fig. 2 Comparison of convergence histories for the Linearized convection model and the 
pressure relaxation procedure for a laminar, incompressible flow in an S-duct. 
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Figure 3. Typical view of the geometry and the grid in the symmetry plane 



Convergence history for an incompressible S-duct 



global sweeps 

Fig 4. Convergence history for a laminar, incompressible flow in an S-duct 
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Figures 5(a-e). Comparison of the laminar axial velocity profiles on the symmetry 

plane at stations 1-5 along the S-duct 
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Figures 6(a-e). Comparison of the computed and measured turbulent axial velocity 

profiles 



Supersonic, Turbulent Convergent-Divergent Nozzle 
- Symmetry plane (x-y) 
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Figure 8. Non-dimensional pressure contours on a symmetry plane 


Figure 9. Velocity vectors on a symmetry plane 


Turbulent, Convergent-Divergent Nozzle 
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Effect of back-pressure (p B ) and divergence angle on the flow 



Figure 11 Comparison of pressure on the lower wall of the symmetry plane for various 

imposed back-pressures 



Effect of back-pressure (p b ) and divergence angle on the flow 
Grid 51x33x33 , Re=500 


Comparison of skin friction coeff C, on the lower wall of the mid-plane 
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Figure 12 Comparison of skin friction coefficient Cj-on the lower wall of the symmetry 

plane for various imposed back-pressures 
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Comparison of axial velocity "u” at the exit (mid-plane) 
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Figure 13 Comparison of axial velocity "u" at the exit cross-plane for various 

imposed back-pressures 







Figure 15 Streamlines on the (x-y) symmetry plane for the imposed back-pressure of 

Pb=6.25p e & Ag/Athroat of 3.51 
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Figure 16 Streamlines on the (x-y) symmetry plane for the imposed back- 
pressure of pb=3.37p e & of 2.3 


Compulation of compressible flow in a C-D nozzle on 



Figure 17 Skin friction coefficient Cf on the lower wall of the (x-y) symmetry plane 
for 31x17x17, 31x33x33 and 51x33x33 grids 
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Computation of compressible flow in a C-D nozzle on 



Figure 20 Cross-flow velocity "v" at the throat for 31x17x17, 31x33x33 and 

51x33x33 grids 
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Figure 21. Typical view of the grid in the symmetry plane 
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Re=5x 1 0 s . M in)eI =3.0 

- Axial velocity "u“ profiles at "station 1" along the symmetry plane 

- for three grids 41x41x41 (Ay min =0.002), 41x41x41 (Ay min =0.0C08) ,61x41x41 (Ay m(n =0.0008) 
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■ Cross flow velocity ‘V profiles at "station 1 " along the symmetry plane 
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Figure 22(a~b). 
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"u " and "v" velocity profiles at ft station 1 " along the mid-plane 
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__ inlet 

- Axial velocity "u" profiles at "station 2" along the symmetry plane 

’ for three grids 41x41x41 (Ay min =0.002), 4 1 x4 1x4 1 (Ay min =0.0008). 61x41x41 (Ay min =0.0008) 
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Turbulent, Supersonic flow in a generic inlet 
Re=5x 10 s . M inlc =3.0 

• Cross flow velocity "v" profiles at ''station 2" along the symmetry plane 
■ for three grids 41x41x41 (Ay min = 0.002), 41x41x41 (Ay min =0.0008), 61x41x41 (Ay n) , n =0.0008) 
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Figure 22(c-d). "u" and "v" velocity profiles at "station 2" along the mid-plane 
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